O termo correlação se refere ao relacionamento entre variáveis e é medida pelo coeficiente de correlação. A correlação simples mede o grau de associação linear entre duas variáveis (\(r_{xy}\)). A correlação parcial mede o grau de associação linear entre duas variáveis mantendo-se outras variáveis constantes. A ordem da correlação parcial vai depender do número de variáveis que ficam constantes. Com uma variável constante (\(r_{xy.z}\)) Z, a correlação parcial de x e y é dita de ordem um. A correlação múltipla mede o grau de associação entre uma variável e um conjunto de variáveis (\(r_{x_1.x_2 x_3 \dots x_p}\)). A correlação canônica mede o grau de associação entre dois grupos de variáveis, com um grupo denotado normalmente por Y e o outro por X. Estatiscamente, a correlação canônica procura testar se os dois conjuntos de variáveis são independentes.
É uma técnica de Análise Multivariada usada para quantificar o grau de correlação (dependência) entre dois grupos de variáveis. Exemplos:
Um exemplo disponível no site da UCLA (https://stats.oarc.ucla.edu/r/dae/canonical-correlation-analysis/) mostra que uma pesquisadora coletou dados sobre três variáveis psicológicas, quatro variáveis acadêmicas (resultados de testes padronizados) e gênero para 600 calouros universitários. Ela está interessada em como o conjunto de variáveis psicológicas se relaciona com as variáveis acadêmicas e de gênero. Em particular, a pesquisadora está interessada em quantas dimensões (variáveis canônicas) são necessárias para entender a associação entre os dois conjuntos de variáveis.
O arquivo de dados (mmreg.csv), possui 600 observações em oito variáveis. As variáveis psicológicas são locus_of_control, auto_conceito e motivação. As variáveis acadêmicas são testes padronizados em leitura (read), escrita (write), matemática (math) e ciências (science). Além disso, a variável female é uma variável binária (zero-um) com 1 para a estudante do sexo feminino.
#Pacotes utilizados
library(tidyverse)
library(skimr)
library(ggplot2)
library(GGally)
library(CCA)
library(CCP)
library(psych)
dados <- read.csv("https://stats.idre.ucla.edu/stat/data/mmreg.csv")
colnames(dados) <- c("Control", "Concept", "Motivation", "Read", "Write", "Math",
"Science", "Gender")
attach(dados)
#Visualizaçao dos dados
glimpse(dados)
## Rows: 600
## Columns: 8
## $ Control <dbl> -0.84, -0.38, 0.89, 0.71, -0.64, 1.11, 0.06, -0.91, 0.45, 0…
## $ Concept <dbl> -0.24, -0.47, 0.59, 0.28, 0.03, 0.90, 0.03, -0.59, 0.03, 0.…
## $ Motivation <dbl> 1.00, 0.67, 0.67, 0.67, 1.00, 0.33, 0.67, 0.67, 1.00, 0.67,…
## $ Read <dbl> 54.8, 62.7, 60.6, 62.7, 41.6, 62.7, 41.6, 44.2, 62.7, 62.7,…
## $ Write <dbl> 64.5, 43.7, 56.7, 56.7, 46.3, 64.5, 39.1, 39.1, 51.5, 64.5,…
## $ Math <dbl> 44.5, 44.7, 70.5, 54.7, 38.4, 61.4, 56.3, 46.3, 54.4, 38.3,…
## $ Science <dbl> 52.6, 52.6, 58.0, 58.0, 36.3, 58.0, 45.0, 36.3, 49.8, 55.8,…
## $ Gender <int> 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1,…
head(dados)
## Control Concept Motivation Read Write Math Science Gender
## 1 -0.84 -0.24 1.00 54.8 64.5 44.5 52.6 1
## 2 -0.38 -0.47 0.67 62.7 43.7 44.7 52.6 1
## 3 0.89 0.59 0.67 60.6 56.7 70.5 58.0 0
## 4 0.71 0.28 0.67 62.7 56.7 54.7 58.0 0
## 5 -0.64 0.03 1.00 41.6 46.3 38.4 36.3 1
## 6 1.11 0.90 0.33 62.7 64.5 61.4 58.0 1
tail(dados)
## Control Concept Motivation Read Write Math Science Gender
## 595 0.46 0.03 1.00 52.1 56.7 62.8 47.1 1
## 596 0.94 -0.30 1.00 60.1 67.1 52.4 55.3 1
## 597 0.23 0.03 1.00 65.4 56.7 65.4 58.0 1
## 598 0.46 0.03 1.00 65.4 51.5 61.4 60.7 1
## 599 0.51 0.03 1.00 54.8 54.1 66.4 41.7 1
## 600 0.25 0.03 0.67 49.5 51.5 55.5 44.4 1
#Estatistica descritiva dos dados
summary(dados)
## Control Concept Motivation Read
## Min. :-2.23000 Min. :-2.620000 Min. :0.0000 Min. :28.3
## 1st Qu.:-0.37250 1st Qu.:-0.300000 1st Qu.:0.3300 1st Qu.:44.2
## Median : 0.21000 Median : 0.030000 Median :0.6700 Median :52.1
## Mean : 0.09653 Mean : 0.004917 Mean :0.6608 Mean :51.9
## 3rd Qu.: 0.51000 3rd Qu.: 0.440000 3rd Qu.:1.0000 3rd Qu.:60.1
## Max. : 1.36000 Max. : 1.190000 Max. :1.0000 Max. :76.0
## Write Math Science Gender
## Min. :25.50 Min. :31.80 Min. :26.00 Min. :0.000
## 1st Qu.:44.30 1st Qu.:44.50 1st Qu.:44.40 1st Qu.:0.000
## Median :54.10 Median :51.30 Median :52.60 Median :1.000
## Mean :52.38 Mean :51.85 Mean :51.76 Mean :0.545
## 3rd Qu.:59.90 3rd Qu.:58.38 3rd Qu.:58.65 3rd Qu.:1.000
## Max. :67.10 Max. :75.50 Max. :74.20 Max. :1.000
skim(dados)
| Name | dados |
| Number of rows | 600 |
| Number of columns | 8 |
| _______________________ | |
| Column type frequency: | |
| numeric | 8 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| Control | 0 | 1 | 0.10 | 0.67 | -2.23 | -0.37 | 0.21 | 0.51 | 1.36 | ▁▂▃▇▃ |
| Concept | 0 | 1 | 0.00 | 0.71 | -2.62 | -0.30 | 0.03 | 0.44 | 1.19 | ▁▁▃▇▃ |
| Motivation | 0 | 1 | 0.66 | 0.34 | 0.00 | 0.33 | 0.67 | 1.00 | 1.00 | ▂▃▁▆▇ |
| Read | 0 | 1 | 51.90 | 10.10 | 28.30 | 44.20 | 52.10 | 60.10 | 76.00 | ▂▇▇▆▂ |
| Write | 0 | 1 | 52.38 | 9.73 | 25.50 | 44.30 | 54.10 | 59.90 | 67.10 | ▁▃▅▆▇ |
| Math | 0 | 1 | 51.85 | 9.41 | 31.80 | 44.50 | 51.30 | 58.37 | 75.50 | ▃▇▇▅▂ |
| Science | 0 | 1 | 51.76 | 9.71 | 26.00 | 44.40 | 52.60 | 58.65 | 74.20 | ▁▅▆▇▂ |
| Gender | 0 | 1 | 0.54 | 0.50 | 0.00 | 0.00 | 1.00 | 1.00 | 1.00 | ▇▁▁▁▇ |
Abaixo está uma lista de alguns métodos de análise que podem ser usadas com estes dados:
Análise de correlação canônica;
Regressões simples separadas – analisar esses dados usando análises de regressão simples separadas para cada variável em um conjunto. As regressões não produzirão resultados multivariados e não relatam informações relativas à dimensionalidade;
A regressão múltipla multivariada.
Os dados para análise de correlações canônicas são formados por dois conjuntos de variáveis, formando dois vetores de variáveis X e Y, que podem ser organizados como demonstrado abaixo.
| OBS | X | Y | ||||||
|---|---|---|---|---|---|---|---|---|
| X1 | X2 | ... | XP | Y1 | Y2 | ... | YP | |
| 1 | ||||||||
| 2 | ||||||||
| 3 | ||||||||
| 4 | ||||||||
| ... | ||||||||
| N |
com estes dados pode-se ter uma matriz de correlaçoes dada por
\[ \mathbf{R}=\left[\begin{array}{ccc} R_{11} & \vdots & R_{12} \\ \dots & \vdots & \dots \\ R_{21} & \vdots & R_{22} \\ \end{array} \right] \]
que pode ser particionada em 4 partes onde \(R_{11}\) é uma matriz de correlações do primeiro grupo; \(R_{22}\) é uma matriz de correlações do segundo grupo; e \(R_{12}\) e \(R_{21}\) são correlações entre as variáveis de um grupo com as do outro grupo, sendo que uma é a transposta da outra.
psych <- dados[, 1:3]
acad <- dados[, 4:8]
ggpairs(psych)
ggpairs(acad)
# Correlações entre os dois conjuntos de variáveis
matcor(psych, acad)
## $Xcor
## Control Concept Motivation
## Control 1.0000000 0.1711878 0.2451323
## Concept 0.1711878 1.0000000 0.2885707
## Motivation 0.2451323 0.2885707 1.0000000
##
## $Ycor
## Read Write Math Science Gender
## Read 1.00000000 0.6285909 0.6792757 0.6906929 -0.04174278
## Write 0.62859089 1.0000000 0.6326664 0.5691498 0.24433183
## Math 0.67927568 0.6326664 1.0000000 0.6495261 -0.04821830
## Science 0.69069291 0.5691498 0.6495261 1.0000000 -0.13818587
## Gender -0.04174278 0.2443318 -0.0482183 -0.1381859 1.00000000
##
## $XYcor
## Control Concept Motivation Read Write Math
## Control 1.0000000 0.17118778 0.24513227 0.37356505 0.35887684 0.3372690
## Concept 0.1711878 1.00000000 0.28857075 0.06065584 0.01944856 0.0535977
## Motivation 0.2451323 0.28857075 1.00000000 0.21060992 0.25424818 0.1950135
## Read 0.3735650 0.06065584 0.21060992 1.00000000 0.62859089 0.6792757
## Write 0.3588768 0.01944856 0.25424818 0.62859089 1.00000000 0.6326664
## Math 0.3372690 0.05359770 0.19501347 0.67927568 0.63266640 1.0000000
## Science 0.3246269 0.06982633 0.11566948 0.69069291 0.56914983 0.6495261
## Gender 0.1134108 -0.12595132 0.09810277 -0.04174278 0.24433183 -0.0482183
## Science Gender
## Control 0.32462694 0.11341075
## Concept 0.06982633 -0.12595132
## Motivation 0.11566948 0.09810277
## Read 0.69069291 -0.04174278
## Write 0.56914983 0.24433183
## Math 0.64952612 -0.04821830
## Science 1.00000000 -0.13818587
## Gender -0.13818587 1.00000000
A matriz \(R_{12}\) mostra as correlações entre as variáveis dos dois grupos. Com p variáveis no vetor X e q variáveis no vetor Y tem-se pq correlações. A análise do relacionamento entre os dois grupos de variáveis por meio das pq correlações é inviável. O objetivo da Análise de Correlação Canônica é resumir o padrão de associação entre os X’s e os Y’s em termos de poucas correlações.
A idéia básica da técnica é:
Calcula-se, inicialmente, 2 combinações lineares, uma de cada conjunto de variáveis, de forma que a correlação entre elas seja máxima. As combinações lineares são denominadas variáveis canônicas e os pares são denominados pares canônicos. O primeiro par deve apresentar a maior correlação entre as variáveis canônicas;
Em seguida, calcula-se 2 outras variáveis canônicas, ou seja, outro par canônico, com a condição de serem ortogonais às primeiras e terem o máximo coeficiente de correlação. O segundo par deve ter a segunda maior correlação entre as variáveis canônicas.
Assim, sucessivamente.
Em geral, o número de dimensões canônicas (pares canônicos) é igual ao número de variáveis no conjunto menor \(min(p,q)\). No entanto, o número de dimensões significativas pode ser ainda menor. O coeficiente de correlação canônica é o coeficiente de correlação de Pearson, em valor absoluto, entre as variáveis canônicas de cada par canônico. Existem, então, \(k=min(p,q)\) pares canônicos e k coeficientes de correlação canônica sendo interpretados somente aqueles estatísticamente significativos.
Sejam \(U_1\) e \(V_1\) o primeiro par canônico, dado por
\[\begin{cases} U_1=a_{11}X_1+a_{21}X_2+a_{31}X_3+ \dots +a_{p1}X_p \\ V_1=b_{11}Y_1+b_{21}Y_2+b_{31}Y_3+ \dots +b_{q1}Y_q \\ \end{cases}\]o problema passa a ser encontrar os vetores \(a_1\) e \(b_1\) que maximizem a correlação entre \(U_1\) e \(V_1\).
Sejam \(U_2\) e \(V_2\) o segundo par canônico, dado por
\[\begin{cases} U_2=a_{12}X_1+a_{22}X_2+a_{32}X_3+ \dots +a_{p2}X_p \\ V_2=b_{12}Y_1+b_{22}Y_2+b_{32}Y_3+ \dots +b_{q2}Y_q \\ \end{cases}\]o problema passa a ser agora encontrar os vetores \(a_2\) e \(b_2\) que maximizem a correlação entre \(U_2\) e \(V_2\) e mantenha a condição de ortogonalidade entre \(U_1\) e \(U_2\) e entre \(V_1\) e \(V_2\).
Considerando as variáveis padronizadas, que é o mesmo que trabalhar com a matriz de correlações, tem-se que combinações lineares das variáveis X são dadas por \(U=a'X\), com \(var(U)=var(a'X)=a'R_{11}a\) e combinações lineares de Y são dadas por \(V=b'Y\), com \(var(V)=var(b'Y)=b'R_{22}b\) e \(cov(U,Y)=cov(a'X,b'Y)=a'R_{12}b\).
Sendo
\[ r_{UV}=\frac{COV(U,V)}{\sqrt{var(U)}\sqrt{var(V)}} \]
pode-se definir \(r_{UV}^2\) como
\[ r_{UV}^2=\frac{[COV(U,V)]^2}{var(U) var(V)}=\frac{(a'R_{12}b)^2}{(a'R_{11}a)(b'R_{22}b)} \]
Assim, o problema pode ser definido como
\[ maximizar \quad (a'R_{12}b) \]
sujeito a
\[ a'R_{11}a=1 \] \[ b'R_{22}b=1 \] \[ a'R_{11}a^*=0 \quad a \neq a^* \] \[ b'R_{22}b^*=0 \quad b \neq b^* \]
Resumidamente, o problema de maximização condicionado encontra os valores de \(a_k\) e \(b_k\) como soluções do seguinte sistema de equações
\[\begin{cases} (R_{12}R_{22}^{-1}R_{21}-\lambda_k R_{11})a_k=0 \\ (R_{21}R_{11}^{-1}R_{12}-\lambda_k R_{22})b_k=0 \\ \end{cases}\]em que \(\lambda_k\) são raízes características e \(a_k\) e \(b_k\) são vetores característicos. Assim, os \(\lambda_k\) satisfazem as seguintes equações características
\[\begin{cases} |R_{12}R_{22}^{-1}R_{21}-\lambda_k R_{11}|=0 \\ |R_{21}R_{11}^{-1}R_{12}-\lambda_k R_{22}|=0 \\ \end{cases}\]ou
\[\begin{cases} |R_{11}^{-1}R_{12}R_{22}^{-1}R_{21}-\lambda_k I|=0 \\ |R_{22}^{-1}R_{21}R_{11}^{-1}R_{12}-\lambda_k I|=0 \\ \end{cases}\]em que, \(\lambda_k\) é a k-ésima maior raiz característica da matriz \(R_{11}^{-1}R_{12}R_{22}^{-1}R_{21}\) ou da matriz \[R_{22}^{-1}R_{21}R_{11}^{-1}R_{12}\].
O vetor “a” é o vetor característico da matriz \(R_{11}^{-1}R_{12}R_{22}^{-1}R_{21}\)
Dado “a”, o vetor “b” pode ser obtido pela relação \(b=\frac{1}{\sqrt{\lambda}}R_{22}^{-1}R_{21}a\).
O vetor “b” pode ser obtido pelas raízes características da matriz \(R_{22}^{-1}R_{21}R_{11}^{-1}R_{12}\).
De forma semelhante, dado “b”, o vetor “a” pode ser calculado por \(a=\frac{1}{\sqrt{\lambda}}R_{11}^{-1}R_{12}b\).
As soluções para as raízes de \(R_{11}^{-1}R_{12}R_{22}^{-1}R_{21}\) ou de \[R_{22}^{-1}R_{21}R_{11}^{-1}R_{12}\] são idênticas.
Tem-se, assim, a solução para os coeficientes das variáveis canônicas que são os elementos dos vetores “a” e “b”. Observe que:
Para cada raiz existe um par de vetores cujos elementos são os coeficientes das variáveis canônicas. Estes coeficientes são chamados pesos canônicos e podem ser calculados na forma bruta ou padronizada.
O coeficiente de correlação canônica é o coeficiente de correlação, em valor absoluto, entre \(U_k\) e \(V_k\). Ele é igual a raiz quadrada da raiz característica, ou seja, \(r_1=\sqrt{\lambda_1}\), \(r_2=\sqrt{\lambda_2}\), etc.
# Calcular e mostrar as correlaçoes canônicas
cc1 <- cc(psych, acad)
cc1$cor #mostra as correlaçoes
## [1] 0.4640861 0.1675092 0.1039911
cc1[3:4] #mostra os coeficientes canônicos brutos
## $xcoef
## [,1] [,2] [,3]
## Control -1.2538339 -0.6214776 -0.6616896
## Concept 0.3513499 -1.1876866 0.8267210
## Motivation -1.2624204 2.0272641 2.0002283
##
## $ycoef
## [,1] [,2] [,3]
## Read -0.044620600 -0.004910024 0.021380576
## Write -0.035877112 0.042071478 0.091307329
## Math -0.023417185 0.004229478 0.009398182
## Science -0.005025152 -0.085162184 -0.109835014
## Gender -0.632119234 1.084642326 -1.794647036
A dimensão 1 teve uma correlação canônica de 0,46 entre os conjuntos de variáveis, enquanto para a dimensão 2 a correlação canônica foi bem menor em 0,17. A correlação canônica da terceira dimensão foi 0,10.
Os coeficientes canônicos brutos são interpretados da mesma forma que os coeficientes da análise de regressão, ou seja, para a variável read, um aumento de uma unidade leva a uma diminuição de 0,0446 na primeira variável canônica do segundo conjunto de dados, tudo o mais constante. Outro exemplo: ser mulher leva a uma diminuição de 0,6321 na dimensão 1 para o conjunto de dados “acadêmico” com tudo mais constante.
# Calcular as cargas canônicas
cc2 <- comput(psych, acad, cc1)
# Mostrar as cargas canônicas
cc2[3:6]
## $corr.X.xscores
## [,1] [,2] [,3]
## Control -0.90404631 -0.3896883 -0.1756227
## Concept -0.02084327 -0.7087386 0.7051632
## Motivation -0.56715106 0.3508882 0.7451289
##
## $corr.Y.xscores
## [,1] [,2] [,3]
## Read -0.3900402 -0.06010654 0.01407661
## Write -0.4067914 0.01086075 0.02647207
## Math -0.3545378 -0.04990916 0.01536585
## Science -0.3055607 -0.11336980 -0.02395489
## Gender -0.1689796 0.12645737 -0.05650916
##
## $corr.X.yscores
## [,1] [,2] [,3]
## Control -0.419555307 -0.06527635 -0.01826320
## Concept -0.009673069 -0.11872021 0.07333073
## Motivation -0.263206910 0.05877699 0.07748681
##
## $corr.Y.yscores
## [,1] [,2] [,3]
## Read -0.8404480 -0.35882541 0.1353635
## Write -0.8765429 0.06483674 0.2545608
## Math -0.7639483 -0.29794884 0.1477611
## Science -0.6584139 -0.67679761 -0.2303551
## Gender -0.3641127 0.75492811 -0.5434036
As correlações acima são entre variáveis observadas e variáveis canônicas que são conhecidas como cargas canônicas. Essas variáveis canônicas são na verdade um tipo de variável latente, ou seja, não observadas. À semelhança da interpretação dos fatores na Análise Fatorial, as variáveis canônicas podem ser identificadas em termos das variáveis com que elas mais se relacionam. Isto pode ser feito por meio dos coeficientes canônicos ou pesos canônicos ou, preferencialmente, através das correlações entre a variável canônica e as variáveis originais. Estas correlações são denominadas cargas canônicas (canonical loadings) ou correlações estruturais. Por exemplo, se \(U_1\) e \(X_1\) são positiva e altamente correlacionadas, \(U_1\) pode ser considerada um reflexo de \(X_1\). Da mesma forma, se \(V_1\) for negativa e altamente correlacionada com \(Y_1\) e positiva e altamente correlacionada com \(Y_3\), então \(V_1\) reflete \(Y_1\) em sentido contrário e \(Y_3\) no mesmo sentido.
São calculados coeficientes de correlação entre as variáveis canônicas e as variáveis originais padronizadas de forma que se tenha os coeficientes de correlação entre as variáveis canônicas U e as variáveis do grupo X; entre as variáveis canônicas U e as variáveis do grupo Y; entre as variáveis canônicas V e as variáveis do grupo Y e entre as variáveis canônicas V e as variáveis do grupo X.
Quando as variáveis do modelo possuem desvios padrão muito diferentes, os coeficientes padronizados permitem comparações mais fáceis entre as variáveis, calculados como demonstrado abaixo.
# coeficientes canonicos padronizados para psych (desvio padrao da diagonal da matriz de var-cov de psych)
s1 <- diag(sqrt(diag(cov(psych))))
s1 %*% cc1$xcoef
## [,1] [,2] [,3]
## [1,] -0.8404196 -0.4165639 -0.4435172
## [2,] 0.2478818 -0.8379278 0.5832620
## [3,] -0.4326685 0.6948029 0.6855370
# coeficientes canonicos padronizados para acad (desvio padrao da diagonal da matriz de var-cov de acad)
s2 <- diag(sqrt(diag(cov(acad))))
s2 %*% cc1$ycoef
## [,1] [,2] [,3]
## [1,] -0.45080116 -0.04960589 0.21600760
## [2,] -0.34895712 0.40920634 0.88809662
## [3,] -0.22046662 0.03981942 0.08848141
## [4,] -0.04877502 -0.82659938 -1.06607828
## [5,] -0.31503962 0.54057096 -0.89442764
Os coeficientes canônicos padronizados são interpretados de maneira semelhante à interpretação dos coeficientes de analise de regressão com variáveis padronizadas. Por exemplo, considere a variável read, um aumento de um desvio padrão em read leva a uma diminuição de 0,45 desvio padrão na pontuação na primeira variável canônica para o segundo conjunto de dados, tudo mais constante.
Para as variáveis psicológicas, a primeira dimensão canônica é mais fortemente influenciada pelo locus de controle (Control) (-0,84) e, para a segunda dimensão, autoconceito (Concept) (-0,84) e motivação (Motivation) (0,69). Para as variáveis acadêmicas mais gênero, a primeira dimensão foi composta por leitura (reading) (-0,45), escrita (writing) (-0,35) e gênero (gender) (-0,32). Para a segunda dimensão, escrita (writing) (0,41), ciência (science) (-0,83) e gênero (gender) (0,54) foram as variáveis dominantes.
Calculados os vetores de coeficientes a e b e os coeficientes de correlação canônica os próximos passos consistem em realizar e interpretar os testes de hipóteses sobre o número de dimensões canônicas. As fórmulas para cálculo podem ser encontradas nas referencias do pacote CCP (https://cran.r-project.org/web/packages/CCP/CCP.pdf)
# testes de dimensões canonicas
rho <- cc1$cor
## Define o numero de observações, numero de variaveis do primeiro grupo e de variáveis do segundo conjunto de variáveis.
n <- dim(psych)[1]
p <- length(psych)
q <- length(acad)
## Calcula os p-values de diferentes testes estatísticos:
p.asym(rho, n, p, q, tstat = "Wilks")
## Wilks' Lambda, using F-approximation (Rao's F):
## stat approx df1 df2 p.value
## 1 to 3: 0.7543611 11.715733 15 1634.653 0.000000000
## 2 to 3: 0.9614300 2.944459 8 1186.000 0.002905057
## 3 to 3: 0.9891858 2.164612 3 594.000 0.091092180
p.asym(rho, n, p, q, tstat = "Hotelling")
## Hotelling-Lawley Trace, using F-approximation:
## stat approx df1 df2 p.value
## 1 to 3: 0.31429738 12.376333 15 1772 0.000000000
## 2 to 3: 0.03980175 2.948647 8 1778 0.002806614
## 3 to 3: 0.01093238 2.167041 3 1784 0.090013176
p.asym(rho, n, p, q, tstat = "Pillai")
## Pillai-Bartlett Trace, using F-approximation:
## stat approx df1 df2 p.value
## 1 to 3: 0.25424936 11.000571 15 1782 0.000000000
## 2 to 3: 0.03887348 2.934093 8 1788 0.002932565
## 3 to 3: 0.01081416 2.163421 3 1794 0.090440474
Conforme mostrado na tabela acima, o primeiro teste das dimensões canônicas testa se todas as três dimensões são significativas (elas são, F = 11,72). O próximo teste verifica se as dimensões 2 e 3 combinadas são significativas (elas são, F = 2,94). Por fim, o último teste analisa se a dimensão 3, por si só, é significativa (não é). Portanto, as dimensões 1 e 2 devem ser significativas, enquanto a dimensão três não é.
Um outro exemplo é de dados de uma amostra de 50 funcionários de uma firma que procura aperfeiçoar testes que podem revelar o potencial para vendas. Foram avaliadas 3 medidas de performance: \(X_1\)= aumento de vendas; \(X_2\)= rendimento das vendas, e \(X_3\)= vendas para novos clientes. Essas medidas foram transformadas para uma escala em que 100 indica performance média. Cada indivíduo foi submetido a 4 testes com o objetivo de medir criatividade, raciocínio mecânico, raciocínio abstrato e habilidade matemática. As notas destes testes formam as variáveis: \(Y_1\)= nota no teste de criatividade; \(Y_2\)= nota no teste de raciocínio mecânico; \(Y_3\)= nota no teste de raciocínio abstrato, e \(Y_4\)= nota no teste de habilidade matemática (Johnson e Wichern, 2007, p. 536). Este exemplo é semelhante ao de Mingoti (2005, p. 146).
# Pacotes
library(tidyverse)
# Direcionado o R para o Diretorio a ser trabalhado
setwd('/Users/jricardofl/Dropbox/tempecon/multivariada')
# Entrada dos dados
dados <- read.csv2("vendedor.csv", header=TRUE, sep=";", dec=".")
dados <- dados[,-1] #retira a primeira coluna
#attach(dados)
#Visualizaçao dos dados
glimpse(dados)
## Rows: 50
## Columns: 7
## $ X1 <dbl> 93.0, 88.8, 95.0, 101.3, 102.0, 95.8, 95.5, 110.8, 102.8, 106.8, 10…
## $ X2 <dbl> 96.0, 91.8, 100.3, 103.8, 107.8, 97.5, 99.5, 122.0, 108.3, 120.5, 1…
## $ X3 <dbl> 97.8, 96.8, 99.0, 106.8, 103.0, 99.3, 99.0, 115.3, 103.8, 102.0, 10…
## $ Y1 <int> 9, 7, 8, 13, 10, 10, 9, 18, 10, 14, 12, 10, 16, 8, 13, 7, 11, 11, 5…
## $ Y2 <int> 12, 10, 12, 14, 15, 14, 12, 20, 17, 18, 17, 18, 17, 10, 10, 9, 12, …
## $ Y3 <int> 9, 10, 9, 12, 12, 11, 9, 15, 13, 11, 12, 8, 11, 11, 8, 5, 11, 11, 1…
## $ Y4 <int> 20, 15, 26, 29, 32, 21, 25, 51, 31, 39, 32, 31, 34, 34, 34, 16, 32,…
head(dados)
## X1 X2 X3 Y1 Y2 Y3 Y4
## 1 93.0 96.0 97.8 9 12 9 20
## 2 88.8 91.8 96.8 7 10 10 15
## 3 95.0 100.3 99.0 8 12 9 26
## 4 101.3 103.8 106.8 13 14 12 29
## 5 102.0 107.8 103.0 10 15 12 32
## 6 95.8 97.5 99.3 10 14 11 21
tail(dados)
## X1 X2 X3 Y1 Y2 Y3 Y4
## 45 94.8 101.8 99.8 7 16 11 24
## 46 103.5 112.0 110.8 18 13 12 37
## 47 89.5 96.0 97.3 7 15 11 14
## 48 84.3 89.8 94.3 8 8 8 9
## 49 104.3 109.5 106.5 14 12 12 36
## 50 106.0 118.5 105.0 12 16 11 39
#Estatistica descritiva dos dados
summary(dados)
## X1 X2 X3 Y1
## Min. : 81.50 Min. : 87.3 Min. : 94.30 Min. : 1.00
## 1st Qu.: 93.55 1st Qu.: 99.5 1st Qu.: 99.08 1st Qu.: 8.25
## Median :100.65 Median :106.2 Median :103.15 Median :10.00
## Mean : 98.84 Mean :106.6 Mean :102.81 Mean :11.22
## 3rd Qu.:105.05 3rd Qu.:114.8 3rd Qu.:106.45 3rd Qu.:14.00
## Max. :110.80 Max. :122.3 Max. :115.30 Max. :18.00
## Y2 Y3 Y4
## Min. : 5.00 Min. : 5.00 Min. : 9.00
## 1st Qu.:12.00 1st Qu.: 9.00 1st Qu.:21.50
## Median :15.00 Median :11.00 Median :31.50
## Mean :14.18 Mean :10.56 Mean :29.76
## 3rd Qu.:17.00 3rd Qu.:12.00 3rd Qu.:37.00
## Max. :20.00 Max. :15.00 Max. :51.00
skim(dados)
| Name | dados |
| Number of rows | 50 |
| Number of columns | 7 |
| _______________________ | |
| Column type frequency: | |
| numeric | 7 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| X1 | 0 | 1 | 98.84 | 7.34 | 81.5 | 93.55 | 100.65 | 105.05 | 110.8 | ▂▅▃▇▆ |
| X2 | 0 | 1 | 106.62 | 10.12 | 87.3 | 99.50 | 106.25 | 114.75 | 122.3 | ▅▆▇▇▇ |
| X3 | 0 | 1 | 102.81 | 4.71 | 94.3 | 99.08 | 103.15 | 106.45 | 115.3 | ▃▅▇▃▁ |
| Y1 | 0 | 1 | 11.22 | 3.95 | 1.0 | 8.25 | 10.00 | 14.00 | 18.0 | ▁▃▇▅▃ |
| Y2 | 0 | 1 | 14.18 | 3.38 | 5.0 | 12.00 | 15.00 | 17.00 | 20.0 | ▁▃▆▇▃ |
| Y3 | 0 | 1 | 10.56 | 2.14 | 5.0 | 9.00 | 11.00 | 12.00 | 15.0 | ▂▃▇▆▁ |
| Y4 | 0 | 1 | 29.76 | 10.54 | 9.0 | 21.50 | 31.50 | 37.00 | 51.0 | ▅▃▇▆▂ |
#Correlações entre os dados
perform <- dados %>%
dplyr::select(X1, X2, X3) %>%
scale() %>% as_tibble
testes <- dados %>%
dplyr::select(Y1, Y2, Y3, Y4) %>%
scale() %>% as_tibble
ggpairs(perform)
ggpairs(testes)
matcor(perform, testes)
## $Xcor
## X1 X2 X3
## X1 1.0000000 0.9260758 0.8840023
## X2 0.9260758 1.0000000 0.8425232
## X3 0.8840023 0.8425232 1.0000000
##
## $Ycor
## Y1 Y2 Y3 Y4
## Y1 1.0000000 0.5907360 0.1469074 0.4126395
## Y2 0.5907360 1.0000000 0.3859502 0.5745533
## Y3 0.1469074 0.3859502 1.0000000 0.5663721
## Y4 0.4126395 0.5745533 0.5663721 1.0000000
##
## $XYcor
## X1 X2 X3 Y1 Y2 Y3 Y4
## X1 1.0000000 0.9260758 0.8840023 0.5720363 0.7080738 0.6744073 0.9273116
## X2 0.9260758 1.0000000 0.8425232 0.5415080 0.7459097 0.4653880 0.9442960
## X3 0.8840023 0.8425232 1.0000000 0.7003630 0.6374712 0.6410886 0.8525682
## Y1 0.5720363 0.5415080 0.7003630 1.0000000 0.5907360 0.1469074 0.4126395
## Y2 0.7080738 0.7459097 0.6374712 0.5907360 1.0000000 0.3859502 0.5745533
## Y3 0.6744073 0.4653880 0.6410886 0.1469074 0.3859502 1.0000000 0.5663721
## Y4 0.9273116 0.9442960 0.8525682 0.4126395 0.5745533 0.5663721 1.0000000
Segundo Mingoti (2005, p.155), “a Análise de Agrupamentos, também conhecida como Análise de Conglomerados, Classificação ou Cluster, tem como objetivo dividir os elementos da amostra, ou população, em grupos de forma que os elementos pertencentes a um mesmo grupo sejam similares entre si com respeito às variáveis (características) que neles foram medidas, e os elementos em grupos diferentes sejam heterogêneos em relação a estas mesmas características”.
Assim, a idéia básica é maximizar a homogeneidade dentro dos grupos ao mesmo tempo em que se maximiza a heterogeneidade entre os grupos. Para formar grupos de objetos tem-se que medir e comparar a semelhança entre eles. A medida de semelhança fica dificultada quando se considera várias variáveis (características).
No processo de agrupamento, primeiro o conjunto de dados será dividido em grupos com elementos parecidos entre si e diferentes dos outros grupos. Depois os grupos encontrados serão analisados até que efetivamente se encontre a existência de padrões entre eles.
Etapas da Análise de Agrupamento - Clustering
No caso univariado é possível que apenas uma análise visual ou gráfica dos dados seja suficiente para decidir sobre os agrupamentos de observaçoes. Contudo, isto não é muito mais difícil quando se tem um conjunto maior de variáveis. Para realizar as técnicas de cluster é necessário definir o conceito de distância.
Escolher um critério de parecença (similaridade ou dissimilaridade);
Formação dos grupos (escolher o algoritmo de agrupamento);
Definição de número de grupos: pode ser a posteriori, como resultado da Análise ou a priori, dado o conhecimento ou conveniência da Análise;
Validação do agrupamento: critério do pesquisador;
Interpretação e Análise: pode-se fazer estatística descritiva, diferença de médias, etc.
Valor numérico que quantifica o grau de semelhança entre um par de objetos. Pode ser de dois tipos: similaridade ou dissimilaridade. Para variáveis quantitativas, quando se buscas agrupar observações as distâncias são as mais usadas. São medidas de dissimilaridade entre objetos. O coeficiente de distância assume valor máximo para objetos totalmente diferentes e valor zero para dois objetos idênticos considerando todas as variáveis.
Em algumas aplicações deseja-se agrupar variáveis ao invés de observações. Neste caso, o coeficiente de correlação de Pearson é mais comumente utilizado. No caso de variáveis binárias, os dados são novamente rearranjados na forma de uma tabela de contigência, com as variáveis formando as categorias.
Os coeficientes de distância mais comuns têm as seguintes propriedades:
Mínimo Zero: Se \(A=B\), então \(D(A,B)=0\);
Positividade: Se \(A \neq B\), então \(D(A,B)>0\);
Simetria: \(D(A,B) = D(B,A)\)
Desigualdade triangular: \(D(A,B)+D(B,C) \geq D(A,C)\)
Em relação às distâncias sabe-se pelo Teorema de Pitágoras que \(c^2=a^2+b^2\), assim, \(c=\sqrt{a^2+b^2}\), que é a Distância Euclidiana entre dois pontos (A e B). Em termos de variáveis, tem-se:
\[ D_{AB}=\sqrt{(X_{2A}-X_{2B})^2+(X_{1B}-X_{1A})^2} \] matricialmente, a Distância Euclidiana é definida por \(D_{AB}=\sqrt{(X_a-X_b)'(X_a-X_b)}\) e possui as seguintes propriedades:
Base geométrica bem definida;
Invariante com relação à transformação de origem;
Invariante com relação à transformação ortogonal;
Não invariante com relação à transformação de escala;
Não invariante com relação à transformação não ortogonal.
A Distância Euclidiana Quadrática é dada por
\[ D_{AB}^2=\sum_{j=1}^p(x_{ja}-x_{jb})^2 \]
A Distância Euclidiana Ponderada é definida por
\[ D_{AB}=\sqrt{(X_a-X_b)'A(X_a-X_b)} \]
com A sendo uma matriz de ponderação (pesos diferentes de acordo com a variância das variáveis). Se \(A=I\), \(D_{AB}\) é a Distância Euclidiana. Se \(A=S^{-1}\), ou seja, inverso da matriz var-cov, tem-se a Distância de Mahalanobis.
Se \(A=diag(S_i^2)^{-1}\), em que \(S_i^2\) é a variância amostral da i-ésima variável aleatória, leva-se em consideração na ponderação apenas as diferenças de variâncias das variáveis.
Se \(A=diag(1/p)\) tem-se a Distância Euclidiana Média;
A Distância de Minkowsky é dada por
\[ D_{AB}=\Big [\sum_i^p w_i | X_{ia} - X_{ib}| ^\lambda \Big]^{\frac{1}{\lambda}} \]
onde \(w_i\) são pesos de ponderação para as variáveis. Se \(\lambda=2\), tem-se a Distância Euclidiana. Esta distância é menos afetada pela presença de valores discrepantes numa amostra do que a distância Euclidiana. Se \(\lambda=1\), tem-se a Distância city-block também conhecida como Manhattan.
Para observações representadas através de variáveis qualitativas, é necessária a criação de variáveis binárias, as quais assumem o valor 1 se uma característica de interesse está presente, e 0, caso contrário. Dessa forma, para um par de observações (i; k) medidos através de p variáveis binárias, considere a seguinte Tabela de Contigência:
| Observação k | ||||
|---|---|---|---|---|
| 1 | 0 | Total | ||
| Observação i | 1 | a | b | a+b |
| 0 | c | d | c+d | |
| Total | a+c | b+d | p=a+b+c+d |
com base nesta tabela é possível listar algumas medidas de similaridade e a idéia básica de cada uma delas.
| \(\frac{a+d}{p}\) | Pesos iguais para pares 1-1 e 0-0 |
| \(\frac{2(a+d)}{2(a+d)+b+c}\) | Peso em dobro para pares 1-1 e 0-0 |
| \(\frac{a+d}{a+d+2(b+c)}\) | Peso em dobro para pares descasados |
| \(\frac{a}{p}\) | Sem pares 0-0 no numerador |
| \(\frac{a}{a+b+c}\) | Sem pares 0-0 no numerador e no denominador |
| \(\frac{2a}{2a+b+c}\) | Sem pares 0 - 0 no numerador e no denominador. |
| Peso em dobro para pares 1 - 1. | |
| \(\frac{a}{a+2(b+c)}\) | Sem pares 0 - 0 no numerador e no denominador. |
| Peso em dobro para pares descasados. | |
| \(\frac{a}{b+c}\) | Razão entre pares 1 - 1 e pares descasados. |
Organiza-se a tabela de contingência e calcula-se o coeficiente de similaridade para cada par de observações. Os coeficientes são dispostos em uma matriz simétrica n x n denominada matriz de similaridade.
Para Mingoti (2005, p.164), “as técnicas de conglomerados ou clusters são frequentemente classificadas em dois tipos: técnicas hierárquicas e não hierárquicas, sendo que as hierárquicas são classificadas em aglomerativas e divisivas”.
Os Métodos hierárquicos aglomerativos partem do princípio de que no início do processo de agrupamento têm-se n conglomerados (cada elemento é um conglomerado). Em cada passo do algoritmo, os elementos vão sendo agrupados, formando novos conglomerados até o ponto em que todos os elementos estão num único grupo.
Em termos de variabilidade, no estágio inicial é a menor possível e no final, a máxima possível, pois todos os elementos estão agrupados em apenas 1 cluster. Se dois elementos aparecem juntos num mesmo grupo em algum estçgio do processo de agrupamento, permanecerão juntos até o final, ou seja, não podem mais ser separados.
Devido a esta propriedade, denominada hierarquia, se pode construiu o Dendograma. Dendograma é um gráfico em forma de árvore no qual a escala vertical indica o nível de similaridade (ou dissimilaridade). No eixo horizontal tem-se os elementos e no vertical a altura correspondente ao nível em que os elementos foram considerados semelhantes. Quanto maior a altura maior a variabilidade, ou seja, maior a heterogeneidade dos grupos a serem unidos.
A outra técnica de construção dos clusters é chamada de Método hierárquivo Divisivo, que se inicia com todos os objetos em um grupo; formam-se subgrupos desagregando os grupos formados até terminar com cada objeto formando um grupo.
Existem ainda os Métodos não hierárquicos, também denominados de métodos de partição. Inicia-se com um número pré definido de grupos e a cada passo procura-se realocar os objetos de maneira a encontrar a melhor partição, isto é, a que minimiza a variância dentro do grupo e maximiza a variância entre grupos. Os métodos não hierárquicos não geram dendrograma, são aplicados a casos (observações) e não a variáveis e são mais indicados para amostras grandes.
\[ D_{I,II} = min \quad ({d_{ij}, i \in I, j \in II}) \]
\[ D_{I,II} = max \quad ({d_{ij}, i \in I, j \in II}) \]
\[ D_{i,j}=\frac{\sum_{i=1}^{n_I} \sum_{j=1}^{n_{II}}d_{ij}}{n_I n_{II}} \]
\[ D_{AB}=(\bar X_a - \bar X_b)'(\bar X_a - \bar X_b) \]
Em um estágio com G grupos, define-se \(SQE\) como a soma \(SQE=SQE_1+SQE_2+SQE_3+ \dots SQE_G\). Quando se junta dois grupos, SQE aumenta. Em cada etapa, a união de todos os pares de grupos é analisada. Junta-se os dois grupos que resultam no menor aumento de \(SQE\)
Segundo Mingoti (2005), “os métodos de ligação simples, completa e da média podem ser utilizados tanto para variáveis quantitativas, quanto qualitativas, ao contrário dos métodos do centróide e de Ward que são apropriados apenas para variáveis quantitativas, já que têm como base a comparação de vetores de médias”.
#Direcionado o R para o Diretorio a ser trabalhado
setwd('/Users/jricardofl/Dropbox/tempecon/dados_gescilene')
#Lendo os dados no R
library(foreign)
library(car)
library(tidyverse)
#library(Rcmdr)
library(corrplot)
library(MVar.pt)
library(cluster)
library(clValid) # Avaliação dos grupos
library(e1071) # Fuzzy K-médias
library(factoextra) # Vizualização de grupos
library(gridExtra) # Ferramentas gráficas
library(ggforce) # Ferramentas gráficas
#Entrada de Dados
dados <- read.dta("dados_gesc.dta")
str(dados)
## 'data.frame': 85 obs. of 16 variables:
## $ x1 : num 60 56 56 61 47 65 50 60 78 68 ...
## $ x2 : num 1 2 2 3 3 3 4 3 2 3 ...
## $ x3 : num 140000 86100 47600 90000 205200 ...
## $ x4 : num 15714 10500 9000 16667 17333 ...
## $ x5 : num 17376 0 0 0 26064 ...
## $ x6 : num 4502 1350 2400 2202 3376 ...
## $ x7 : num 26400 21600 20400 18000 22800 21600 24000 13200 16800 21600 ...
## $ x8 : num 48878 25800 26080 22902 55204 ...
## $ x9 : num 2 0 0 0 3 0 2 0 2 0 ...
## $ x10: num 0.111 0 0 0.222 0.111 ...
## $ x11: num 0.143 0 0.143 0 0 ...
## $ x12: num 20 0 1 5 20 20 10 15 20 10 ...
## $ x13: num 2 0 1 2 2 2 2 1 1 1 ...
## $ x14: num 1 0 1 0 0 0 0 0 0 0 ...
## $ x15: num 0.26 0.12 0.09 0.06 0.06 ...
## $ x16: num 1 0 1 0 1 0 1 1 1 1 ...
## - attr(*, "datalabel")= chr ""
## - attr(*, "time.stamp")= chr ""
## - attr(*, "formats")= chr [1:16] "%9.0g" "%9.0g" "%9.0g" "%9.0g" ...
## - attr(*, "types")= int [1:16] 254 254 254 254 254 254 254 254 254 254 ...
## - attr(*, "val.labels")= chr [1:16] "" "" "" "" ...
## - attr(*, "var.labels")= chr [1:16] "" "" "" "" ...
## - attr(*, "version")= int -7
# Nomes das Variáveis utilizadas
#x1 Idade
#x2 Escolaridade
#x3 Renda Bruta Total
#x4 Produtividade
#x5 Valor total anual de mão de obra permanente
#x6 Valor Total com custo anual com insumos agrícolas nas ativ. Irrig. em 2014
#x7 Custo com energia eletrica e agua
#x8 Capital total empregado na ativ irrig em 2014
#x9 N total de empregados
#x10 Indice de introducao de inovacao (III)
#x11 Indice de Inovacoes realizadas (II)
#x12 Gastos para desenvolver as atividade de inovacao sobre faturamento de 2014
#x13 Quantidade de tecnologia agricola utilizada na atividade irrigada
#x14 Dummy se a empresa realizou atividades de treinamento e capacitacao de recursos humanos entre 2010 a 2014
#x15 Indice de Fonte de Informacao
#x16 Dummy se empresa esteve envolvida em atividades cooperativas entre os anos de 2010 e 2014
#####################################################
#dados <- dados1 %>% select(-x14, -x16)
dados <- dados[,-c(14,16)]
#Se quiser padronizar os dados
dados.pad <-as.data.frame(scale(dados))
#Analise Fatorial
fatorial1 <- factanal(dados.pad, factors=4, rotation="none", na.action=na.omit)
fatorial2 <- factanal(dados.pad, factors=4, rotation="varimax", na.action=na.omit, scores = c("regression"))
#Analise de Cluster
# Função que obtém o WSS para o método hierárquico
# NÃO MUDAR!
get_wss <- function(d, cluster){
d <- stats::as.dist(d)
cn <- max(cluster)
clusterf <- as.factor(cluster)
clusterl <- levels(clusterf)
cnn <- length(clusterl)
if (cn != cnn) {
warning("cluster renumbered because maximum != number of clusters")
for (i in 1:cnn) cluster[clusterf == clusterl[i]] <- i
cn <- cnn
}
cwn <- cn
# Compute total within sum of square
dmat <- as.matrix(d)
within.cluster.ss <- 0
for (i in 1:cn) {
cluster.size <- sum(cluster == i)
di <- as.dist(dmat[cluster == i, cluster == i])
within.cluster.ss <- within.cluster.ss + sum(di^2)/cluster.size
}
within.cluster.ss
}
# Função que constrói o "gráfico do cotovelo" e aponta o K ótimo
# NÃO MUDAR!
elbow.plot <- function(x, kmax = 15, alg = "kmeans") {
# alg = c("kmeans", "cmeans", "hclust")
wss <- c()
if (alg == "kmeans") {
for (i in 1:kmax) {
set.seed(13)
tmp <- kmeans(x, i)
# wss[i] <- get_wss(dist(x), tmp$cluster)
wss[i] <- tmp$tot.withinss
}
tmp <- data.frame(k = 1:kmax, wss)
max_k <- max(tmp$k)
max_k_wss <- tmp$wss[which.max(tmp$k)]
max_wss <- max(tmp$wss)
max_wss_k <- tmp$k[which.max(tmp$wss)]
max_df <- data.frame(x = c(max_wss_k, max_k), y = c(max_wss, max_k_wss))
tmp_lm <- lm(max_df$y ~ max_df$x)
d <- c()
for(i in 1:kmax) {
d <- c(d, abs(coef(tmp_lm)[2]*i - tmp$wss[i] + coef(tmp_lm)[1]) /
sqrt(coef(tmp_lm)[2]^2 + 1^2))
}
tmp$d <- d
ggplot(data = tmp, aes(k, wss)) +
geom_line() +
geom_segment(aes(x = k[1], y = wss[1],
xend = max(k), yend = wss[which.max(k)]),
linetype = "dashed") +
geom_point(aes(size = (d == max(d)), color = (d == max(d))),
show.legend = FALSE) +
scale_size_manual(values = c(2,5)) +
scale_color_manual(values = c("black", "red")) +
labs(x = "Number of clusters",
y = "Total within-cluster sum of squares",
title = "Elbow plot for the K-means method") +
theme_bw()
}
else if (alg == "cmeans") {
for (i in 1:kmax) {
if (i == 1) {
wss[i] <- get_wss(dist(x), rep(1, nrow(x)))
}
else {
set.seed(13)
tmp <- cmeans(x, i)
wss[i] <- get_wss(dist(x), tmp$cluster)
# wss[i] <- tmp$sumsqrs$tot.within.ss
}
}
tmp <- data.frame(k = 1:kmax, wss)
max_k <- max(tmp$k)
max_k_wss <- tmp$wss[which.max(tmp$k)]
max_wss <- max(tmp$wss)
max_wss_k <- tmp$k[which.max(tmp$wss)]
max_df <- data.frame(x = c(max_wss_k, max_k), y = c(max_wss, max_k_wss))
tmp_lm <- lm(max_df$y ~ max_df$x)
d <- c()
for(i in 1:kmax) {
d <- c(d, abs(coef(tmp_lm)[2]*i - tmp$wss[i] + coef(tmp_lm)[1]) /
sqrt(coef(tmp_lm)[2]^2 + 1^2))
}
tmp$d <- d
ggplot(data = tmp, aes(k, wss)) +
geom_line() +
geom_segment(aes(x = k[1], y = wss[1],
xend = max(k), yend = wss[which.max(k)]),
linetype = "dashed") +
geom_point(aes(size = (d == max(d)), color = (d == max(d))),
show.legend = FALSE) +
scale_size_manual(values = c(2,5)) +
scale_color_manual(values = c("black", "red")) +
labs(x = "Number of clusters",
y = "Total within-cluster sum of squares",
title = "Elbow plot for the fuzzy K-means method") +
theme_bw()
}
else if (alg == "hclust") {
for (i in 1:kmax) {
set.seed(13)
tmp <- hcut(x, i)
wss[i] <- get_wss(dist(x), tmp$cluster)
}
tmp <- data.frame(k = 1:kmax, wss)
max_k <- max(tmp$k)
max_k_wss <- tmp$wss[which.max(tmp$k)]
max_wss <- max(tmp$wss)
max_wss_k <- tmp$k[which.max(tmp$wss)]
max_df <- data.frame(x = c(max_wss_k, max_k), y = c(max_wss, max_k_wss))
tmp_lm <- lm(max_df$y ~ max_df$x)
d <- c()
for(i in 1:kmax) {
d <- c(d, abs(coef(tmp_lm)[2]*i - tmp$wss[i] + coef(tmp_lm)[1]) /
sqrt(coef(tmp_lm)[2]^2 + 1^2))
}
tmp$d <- d
ggplot(data = tmp, aes(k, wss)) +
geom_line() +
geom_segment(aes(x = k[1], y = wss[1],
xend = max(k), yend = wss[which.max(k)]),
linetype = "dashed") +
geom_point(aes(size = (d == max(d)), color = (d == max(d))),
show.legend = FALSE) +
scale_size_manual(values = c(2,5)) +
scale_color_manual(values = c("black", "red")) +
labs(x = "Number of clusters",
y = "Total within-cluster sum of squares",
title = "Elbow plot for the hierarchical method") +
theme_bw()
}
}
# Função vizualização dos grupos
# NÃO MUDAR!
cluster_viz <- function(data, clusters,
axes = c(1, 2), geom = c("point", "text"), repel = TRUE,
show.clust.cent = TRUE, ellipse = TRUE, ellipse.type = "convex",
ellipse.level = 0.95, ellipse.alpha = 0.2, shape = NULL,
pointsize = 1.5, labelsize = 12, main = "Cluster plot",
ggtheme = theme_bw()) {
require(factoextra)
data <- scale(data)
pca <- stats::prcomp(data, scale = FALSE, center = FALSE)
ind <- facto_summarize(pca, element = "ind", result = "coord", axes = axes)
eig <- get_eigenvalue(pca)[axes, 2]
xlab <- paste0("Dim", axes[1], " (", round(eig[1], 1), "%)")
ylab <- paste0("Dim", axes[2], " (", round(eig[2], 1), "%)")
colnames(ind)[2:3] <- c("x", "y")
label_coord <- ind
lab <- NULL
if ("text" %in% geom)
lab <- "name"
if (is.null(shape))
shape <- "cluster"
plot.data <- cbind.data.frame(ind, cluster = clusters, stringsAsFactors = TRUE)
label_coord <- cbind.data.frame(label_coord, cluster = clusters, stringsAsFactors = TRUE)
p <- ggpubr::ggscatter(plot.data, "x", "y", color = "cluster",
shape = shape, size = pointsize, point = "point" %in% geom,
label = lab, font.label = labelsize, repel = repel,
mean.point = show.clust.cent, ellipse = ellipse, ellipse.type = ellipse.type,
ellipse.alpha = ellipse.alpha, ellipse.level = ellipse.level,
main = main, xlab = xlab, ylab = ylab, ggtheme = ggtheme)
p
}
# Encontrar número ótimo de grupos para o método hierárquico
set.seed(123456789)
p3 <- elbow.plot(dados.pad, alg = "hclust")
p3
# K "ótimo" é igual a 4. Mudar de acordo com o seu resultado
d <-dist(dados.pad, method = "euclidean")
h.fit <-hclust(d, method = "ward.D")
plot(h.fit, main='Dendograma - Método Hierárquico', xlab='Cluster das Observações - Distância Euclidiana e Método de Ward', ylab='Altura')
groups_a <- cutree(h.fit, k=3)
groups_a
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2
## 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
## 2 2 1 2 2 2 2 2 2 1 2 2 2 2 2 2 2 3 2 3 2 2 3 2 2 3
## 79 80 81 82 83 84 85
## 3 2 2 2 2 2 2
groups_b <- cutree(h.fit, k=4)
groups_b
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 3
## 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
## 3 3 1 2 2 2 2 2 3 1 2 2 3 3 2 2 2 4 2 4 2 2 4 3 2 4
## 79 80 81 82 83 84 85
## 4 2 2 2 2 2 2
groups_c <- cutree(h.fit, k=5)
dendogram3 <- rect.hclust(h.fit, k=5, border="red")
dendogram3
## [[1]]
## 2 3 4 6 9 10 11 13 17 18 28 29 31 33 37 42 43 48 49 55 62
## 2 3 4 6 9 10 11 13 17 18 28 29 31 33 37 42 43 48 49 55 62
##
## [[2]]
## 1 5 7 8 12 14 15 16 19 20 21 22 23 24 25 26 27 30 32 34 35 36 38 39 40 41
## 1 5 7 8 12 14 15 16 19 20 21 22 23 24 25 26 27 30 32 34 35 36 38 39 40 41
## 44 45 46 47 51
## 44 45 46 47 51
##
## [[3]]
## 70 72 75 78 79
## 70 72 75 78 79
##
## [[4]]
## 52 53 54 61 65 66 76
## 52 53 54 61 65 66 76
##
## [[5]]
## 50 56 57 58 59 60 63 64 67 68 69 71 73 74 77 80 81 82 83 84 85
## 50 56 57 58 59 60 63 64 67 68 69 71 73 74 77 80 81 82 83 84 85
summary(fatorial2$scores)
## Factor1 Factor2 Factor3 Factor4
## Min. :-1.1249 Min. :-1.11082 Min. :-3.1932 Min. :-2.11638
## 1st Qu.:-0.6559 1st Qu.:-0.21990 1st Qu.:-0.2904 1st Qu.:-0.18303
## Median :-0.3633 Median :-0.11821 Median :-0.1839 Median :-0.14019
## Mean : 0.0000 Mean : 0.00000 Mean : 0.0000 Mean : 0.00000
## 3rd Qu.: 0.3353 3rd Qu.:-0.04689 3rd Qu.: 0.0107 3rd Qu.:-0.06075
## Max. : 3.3937 Max. : 7.06890 Max. : 6.5278 Max. : 7.41904
attach(dados)
by(dados, groups_c, summary)
## groups_c: 1
## x1 x2 x3 x4
## Min. :29.0 Min. :1.000 Min. : 33000 Min. : 6667
## 1st Qu.:49.0 1st Qu.:2.000 1st Qu.: 87150 1st Qu.:13345
## Median :52.0 Median :3.000 Median :127400 Median :16667
## Mean :53.1 Mean :2.903 Mean :147082 Mean :16839
## 3rd Qu.:60.0 3rd Qu.:3.000 3rd Qu.:166400 3rd Qu.:21027
## Max. :70.0 Max. :5.000 Max. :461700 Max. :30000
## x5 x6 x7 x8
## Min. : 0 Min. : 1071 Min. :12000 Min. : 16660
## 1st Qu.: 0 1st Qu.: 1401 1st Qu.:18000 1st Qu.: 30268
## Median :17376 Median : 1900 Median :19200 Median : 43008
## Mean :18577 Mean : 3339 Mean :20923 Mean : 44942
## 3rd Qu.:26064 3rd Qu.: 3378 3rd Qu.:24000 3rd Qu.: 55152
## Max. :69504 Max. :27480 Max. :38400 Max. :116484
## x9 x10 x11 x12
## Min. :0.000 Min. :0.0000 Min. :0.00000 Min. :10.0
## 1st Qu.:0.000 1st Qu.:0.1111 1st Qu.:0.00000 1st Qu.:20.0
## Median :2.000 Median :0.1111 Median :0.00000 Median :20.0
## Mean :2.065 Mean :0.1254 Mean :0.03226 Mean :22.9
## 3rd Qu.:3.000 3rd Qu.:0.1111 3rd Qu.:0.00000 3rd Qu.:25.0
## Max. :8.000 Max. :0.2222 Max. :0.14286 Max. :50.0
## x13 x15
## Min. :0.000 Min. :0.03
## 1st Qu.:1.000 1st Qu.:0.06
## Median :1.000 Median :0.12
## Mean :1.323 Mean :0.12
## 3rd Qu.:2.000 3rd Qu.:0.17
## Max. :2.000 Max. :0.30
## ------------------------------------------------------------
## groups_c: 2
## x1 x2 x3 x4
## Min. :50.00 Min. :1.000 Min. : 21250 Min. : 3333
## 1st Qu.:61.00 1st Qu.:2.000 1st Qu.: 49200 1st Qu.:13333
## Median :65.00 Median :2.000 Median : 86100 Median :16429
## Mean :65.95 Mean :2.095 Mean : 91014 Mean :17585
## 3rd Qu.:71.00 3rd Qu.:2.000 3rd Qu.:126000 3rd Qu.:21000
## Max. :78.00 Max. :3.000 Max. :160402 Max. :31000
## x5 x6 x7 x8
## Min. : 0 Min. : 880 Min. : 9240 Min. :12870
## 1st Qu.: 0 1st Qu.: 1650 1st Qu.:16800 1st Qu.:23410
## Median : 0 Median : 2003 Median :20400 Median :26080
## Mean : 8834 Mean : 3469 Mean :19069 Mean :33256
## 3rd Qu.:17376 3rd Qu.: 3000 3rd Qu.:21600 3rd Qu.:38000
## Max. :34752 Max. :25300 Max. :26400 Max. :62740
## x9 x10 x11 x12
## Min. :0.0000 Min. :0.0000 Min. :0.00000 Min. : 0.000
## 1st Qu.:0.0000 1st Qu.:0.1111 1st Qu.:0.00000 1st Qu.: 5.000
## Median :0.0000 Median :0.1111 Median :0.00000 Median :10.000
## Mean :0.9524 Mean :0.1058 Mean :0.03401 Mean : 8.619
## 3rd Qu.:2.0000 3rd Qu.:0.1111 3rd Qu.:0.00000 3rd Qu.:10.000
## Max. :4.0000 Max. :0.2222 Max. :0.14286 Max. :20.000
## x13 x15
## Min. :0.000 Min. :0.060
## 1st Qu.:1.000 1st Qu.:0.120
## Median :1.000 Median :0.120
## Mean :1.286 Mean :0.149
## 3rd Qu.:2.000 3rd Qu.:0.160
## Max. :2.000 Max. :0.460
## ------------------------------------------------------------
## groups_c: 3
## x1 x2 x3 x4
## Min. :26.00 Min. :1.000 Min. : 105000 Min. :11284
## 1st Qu.:35.00 1st Qu.:5.000 1st Qu.: 208000 1st Qu.:22381
## Median :42.00 Median :5.000 Median : 403580 Median :30000
## Mean :43.24 Mean :4.952 Mean : 567728 Mean :30888
## 3rd Qu.:51.00 3rd Qu.:6.000 3rd Qu.: 720000 3rd Qu.:34714
## Max. :65.00 Max. :8.000 Max. :2380000 Max. :62963
## x5 x6 x7 x8
## Min. : 10080 Min. : 1350 Min. : 3720 Min. : 51280
## 1st Qu.: 42000 1st Qu.: 55000 1st Qu.:12360 1st Qu.:148600
## Median : 83040 Median :145000 Median :15840 Median :309824
## Mean :117334 Mean :146283 Mean :20742 Mean :306951
## 3rd Qu.:167202 3rd Qu.:220000 3rd Qu.:25800 3rd Qu.:420600
## Max. :312000 Max. :350000 Max. :51600 Max. :608400
## x9 x10 x11 x12
## Min. : 0.000 Min. :0.0000 Min. :0.0000 Min. : 0.000
## 1st Qu.: 4.000 1st Qu.:0.2222 1st Qu.:0.0000 1st Qu.: 0.000
## Median : 8.000 Median :0.3333 Median :0.1429 Median : 5.000
## Mean : 9.238 Mean :0.2910 Mean :0.1463 Mean : 5.048
## 3rd Qu.:13.000 3rd Qu.:0.4444 3rd Qu.:0.2857 3rd Qu.:10.000
## Max. :25.000 Max. :0.5556 Max. :0.4286 Max. :10.000
## x13 x15
## Min. :1.000 Min. :0.0000
## 1st Qu.:3.000 1st Qu.:0.4700
## Median :4.000 Median :0.5900
## Mean :3.571 Mean :0.5671
## 3rd Qu.:4.000 3rd Qu.:0.7100
## Max. :6.000 Max. :0.8900
## ------------------------------------------------------------
## groups_c: 4
## x1 x2 x3 x4
## Min. :39.00 Min. :3.000 Min. : 176000 Min. : 8000
## 1st Qu.:49.50 1st Qu.:6.000 1st Qu.: 714000 1st Qu.:13583
## Median :52.00 Median :7.000 Median :1015000 Median :19444
## Mean :52.43 Mean :6.143 Mean :1017429 Mean :19537
## 3rd Qu.:57.00 3rd Qu.:7.000 3rd Qu.:1446500 3rd Qu.:20800
## Max. :63.00 Max. :7.000 Max. :1610000 Max. :40545
## x5 x6 x7 x8
## Min. : 24000 Min. : 49000 Min. :15240 Min. :101280
## 1st Qu.:125544 1st Qu.:134150 1st Qu.:28380 1st Qu.:472710
## Median :242448 Median :235000 Median :37200 Median :605008
## Mean :223483 Mean :244614 Mean :44286 Mean :588103
## 3rd Qu.:279696 3rd Qu.:330000 3rd Qu.:59400 3rd Qu.:739880
## Max. :487452 Max. :500000 Max. :82000 Max. :985252
## x9 x10 x11 x12
## Min. : 0.00 Min. :0.5556 Min. :0.4286 Min. : 1.00
## 1st Qu.: 8.50 1st Qu.:0.5556 1st Qu.:0.5000 1st Qu.: 6.50
## Median :22.00 Median :0.5556 Median :0.6429 Median :10.00
## Mean :20.14 Mean :0.6508 Mean :0.5918 Mean :11.71
## 3rd Qu.:28.50 3rd Qu.:0.7222 3rd Qu.:0.6786 3rd Qu.:15.00
## Max. :45.00 Max. :0.8889 Max. :0.7143 Max. :28.00
## x13 x15
## Min. : 5.000 Min. :0.6700
## 1st Qu.: 6.000 1st Qu.:0.7800
## Median : 6.000 Median :0.8800
## Mean : 7.286 Mean :0.8629
## 3rd Qu.: 7.500 3rd Qu.:0.9650
## Max. :13.000 Max. :1.0000
## ------------------------------------------------------------
## groups_c: 5
## x1 x2 x3 x4
## Min. :27.0 Min. :7.0 Min. :4290000 Min. :18033
## 1st Qu.:29.0 1st Qu.:7.0 1st Qu.:4375000 1st Qu.:25000
## Median :29.0 Median :7.0 Median :6160000 Median :34310
## Mean :35.4 Mean :7.4 Mean :5616500 Mean :29814
## 3rd Qu.:40.0 3rd Qu.:8.0 3rd Qu.:6237500 3rd Qu.:35227
## Max. :52.0 Max. :8.0 Max. :7020000 Max. :36500
## x5 x6 x7 x8
## Min. : 391680 Min. : 553000 Min. : 0 Min. :1829004
## 1st Qu.: 456840 1st Qu.:1200000 1st Qu.: 60000 1st Qu.:1855840
## Median : 920004 Median :1300000 Median : 60000 Median :2747440
## Mean : 809705 Mean :1730600 Mean :105600 Mean :2991657
## 3rd Qu.:1080000 3rd Qu.:1900000 3rd Qu.:144000 3rd Qu.:3376000
## Max. :1200000 Max. :3700000 Max. :264000 Max. :5150000
## x9 x10 x11 x12 x13
## Min. : 0.0 Min. :0.4444 Min. :0.2143 Min. : 0 Min. : 3.0
## 1st Qu.:30.0 1st Qu.:0.7778 1st Qu.:0.4286 1st Qu.: 5 1st Qu.: 6.0
## Median :44.0 Median :0.7778 Median :0.7143 Median : 5 Median : 8.0
## Mean :50.8 Mean :0.7556 Mean :0.6143 Mean : 8 Mean : 7.4
## 3rd Qu.:90.0 3rd Qu.:0.7778 3rd Qu.:0.8571 3rd Qu.:10 3rd Qu.: 8.0
## Max. :90.0 Max. :1.0000 Max. :0.8571 Max. :20 Max. :12.0
## x15
## Min. :0.460
## 1st Qu.:0.760
## Median :0.810
## Mean :0.762
## 3rd Qu.:0.860
## Max. :0.920
# Agrupamento usando o método hierárquico - 2 forma
nclust = 5 # Número de grupos de acordo com os gráficos
fit.hclust <- hcut(dados.pad, 5) # Ajustar K de acordo com o gráfico
fit.hclust
##
## Call:
## stats::hclust(d = x, method = hc_method)
##
## Cluster method : ward.D2
## Distance : euclidean
## Number of objects: 85
g_cluster <- cluster_viz(dados.pad, as.factor(fit.hclust$cluster), geom = "point", main = "Gráfico Cluster para o Método Hierárquico")
plot(g_cluster)
São métodos de partição de “n” objetos em “k” grupos. Os “k” grupos são pré-determinados. A cada passo verifica-se se os objetos estão alocados da “melhor” maneira. Função objetivo: procura minimizar a variância dentro do grupo e maximizar a variância entre os grupos.
Não fornece dendograma. Isto porque se em algum passo do algoritmo dois elementos tiverem sido colocados num mesmo cluster, não necessariamente “estarão juntos” na partição final. Bom quando se tem um grande número de observações. O problema principal aqui não é o número de grupos e sim a melhor forma de alocar os elementos em k grupos. O método mais usado é o das k-médias (k-means).
Resumidamente, o método das k-means é composto de quatro passos. No primeiro são selecionados k centróides para inicializar o processo de partição. Cada elemento do conjunto de dados é, então, comparado com cada centróide inicial, através de uma medida de distância que, em geral, é a Euclidiana. O elemento é alocado ao grupo cuja distância é a menor. No terceiro passo, recalcula-se os valores dos centróides para cada novo grupo formado. Com estes novos centróides, repete-se o passo 2. Isto deve ser repetido até que nenhuma realocação de elementos seja necessária.
cluster_kmeans <- elbow.plot(dados.pad)
cluster_kmeans
kmeans1 <- kmeans(dados.pad, centers=5)
#Centróides
kmeans1$centers
## x1 x2 x3 x4 x5 x6
## 1 -1.34515006 1.8072111 3.701528552 0.8027013 3.2778778 3.23549118
## 2 -0.02424868 1.1862676 0.287208385 -0.1871013 0.5550533 0.17399627
## 3 -0.02004540 -0.3831231 -0.355117495 -0.3418032 -0.3879402 -0.32287084
## 4 -0.79009568 0.6492143 0.007068955 1.1307878 0.1277874 0.01285434
## 5 0.83455884 -0.6890467 -0.397959050 -0.5056032 -0.4355642 -0.31782898
## x7 x8 x9 x10 x11 x12 x13
## 1 2.4501995 3.49247328 2.7059130 2.2639164 2.1490201 -0.5026098 1.8761180
## 2 0.5309120 0.34475417 0.7714212 1.8022454 2.0470493 -0.1411494 1.8303262
## 3 -0.1966927 -0.36414512 -0.3560014 -0.5085221 -0.5070145 1.0175370 -0.5500592
## 4 -0.1810064 0.02647696 0.1454207 0.2944307 0.1338101 -0.7891520 0.3802517
## 5 -0.2694298 -0.37661625 -0.4316555 -0.5572274 -0.4915234 -0.4539516 -0.5032957
## x15
## 1 1.4079157
## 2 1.7417060
## 3 -0.6848525
## 4 0.9214138
## 5 -0.6137018
#Plot dos Clusters
clusplot(dados.pad, kmeans1$cluster, main="2D Representação da solução dos Clusters",
color=TRUE, shade=TRUE, labels=2, lines=0)
# Agrupamento usando o K-means - 2 forma
fit.kmeans <- kmeans(dados.pad, 5) # Ajustar K de acordo com o gráfico
fit.kmeans
## K-means clustering with 5 clusters of sizes 18, 5, 29, 26, 7
##
## Cluster means:
## x1 x2 x3 x4 x5 x6
## 1 -0.79009568 0.6492143 0.007068955 1.1307878 0.1277874 0.01285434
## 2 -1.34515006 1.8072111 3.701528552 0.8027013 3.2778778 3.23549118
## 3 -0.02004540 -0.3831231 -0.355117495 -0.3418032 -0.3879402 -0.32287084
## 4 0.83455884 -0.6890467 -0.397959050 -0.5056032 -0.4355642 -0.31782898
## 5 -0.02424868 1.1862676 0.287208385 -0.1871013 0.5550533 0.17399627
## x7 x8 x9 x10 x11 x12 x13
## 1 -0.1810064 0.02647696 0.1454207 0.2944307 0.1338101 -0.7891520 0.3802517
## 2 2.4501995 3.49247328 2.7059130 2.2639164 2.1490201 -0.5026098 1.8761180
## 3 -0.1966927 -0.36414512 -0.3560014 -0.5085221 -0.5070145 1.0175370 -0.5500592
## 4 -0.2694298 -0.37661625 -0.4316555 -0.5572274 -0.4915234 -0.4539516 -0.5032957
## 5 0.5309120 0.34475417 0.7714212 1.8022454 2.0470493 -0.1411494 1.8303262
## x15
## 1 0.9214138
## 2 1.4079157
## 3 -0.6848525
## 4 -0.6137018
## 5 1.7417060
##
## Clustering vector:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
## 3 4 4 4 3 3 4 4 4 4 4 3 4 3 3 3 4 4 3 3 3 3 3 3 3 3
## 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
## 3 4 4 4 4 3 4 3 3 3 4 3 3 3 4 4 4 3 3 3 3 4 4 3 3 5
## 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
## 5 5 4 1 4 1 1 1 5 4 1 1 5 5 1 1 1 2 1 2 1 1 2 5 1 2
## 79 80 81 82 83 84 85
## 2 1 4 1 1 1 1
##
## Within cluster sum of squares by cluster:
## [1] 86.28616 143.25082 64.37370 48.86236 43.95670
## (between_SS / total_SS = 67.1 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"
g_kmeans <- cluster_viz(dados.pad, as.factor(fit.kmeans$cluster), geom = "point",
main = "Gráfico Cluster para o Método K-médias")
plot(g_kmeans)
A Análise Discriminante Linear (LDA) é uma técnica de Análise Multivariada usada com o objetivo de discriminar (separar, diferenciar) populações e/ou classificar observações em populações pré definidas. Para a sua aplicação, é necessário que os grupos para os quais cada elemento amostral pode ser classificado sejam pré-definidos, ou seja, conhecidos a priori considerando-se suas características gerais. Este conhecimento permite a elaboração de uma função matemática chamada de regra de classificação ou discriminação, que é utilizada para classificar novos elementos amostrais nos grupos já existentes. A definição dos grupos é feita de acordo com o problema de pesquisa e os objetivos do estudo.
A LDA, originalmente desenvolvida por R A Fisher (1936) para classificar indivíduos em um dos dois grupos claramente definidos. Posteriormente, foi expandido para classificar os assuntos em mais de dois grupos. A LDA é uma técnica para reduzir o número de dimensões (ou seja, variáveis) em um conjunto de dados, mantendo o máximo de informações possível. Basicamente, ajuda a encontrar a combinação linear das variáveis originais que fornecem a melhor separação possível entre os grupos.
O objetivo é construir uma regra de classificação, baseada na teoria das probabilidades, que minimize o número de classificações incorretas, ou seja, o erro de dizer que um elemento pertence à população A quando na verdade ele pertence à população B. Uma informação importante é saber quais variáveis mais ajudam a explicar a diferença entre os grupos. A idéia básica é saber se as respectivas médias dos grupos provêm de distintas populações.
A Análise Discriminante tem, assim, dois objetivos básicos:
Discriminar – consiste em encontrar funções capazes de explicar as diferenças entre os grupos com base em várias variáveis; estas funções captam o máximo de separação entre os grupos.
Classificar – consiste no uso das funções para alocar novas observações em um dos grupos pré-estabelecidos. Assim, pode-se dizer que existem dois tipos de Análise Discriminante: Descritiva e Preditiva.
As principais aplicações ou exemplos são:
Prevendo o sucesso ou fracasso de novos produtos;
Aceitar ou rejeitar a admissão de um candidato;
Prevendo a categoria de risco de crédito para uma pessoa;
Classificando os pacientes em diferentes categorias.
Para aplicação da LDA deve-se dispor de uma amostra de cada uma das g populações e dados de p variáveis observadas para cada objeto das amostras que podem ser de tamanhos diferentes. Ao contrário da Análise de Cluster em que se procura formar grupos homogêneos na amostra, na LDA parte-se de um número pré-definido de grupos considerados amostras de diferentes populações.
A base de dados utilizada para o exemplo é o da flor Iris ou conjunto de dados Iris de Fisher, um conjunto de dados multivariados introduzido pelo estatístico e biólogo britânico Ronald Fisher em seu artigo de 1936, “O uso de múltiplas medições em problemas taxonômicos, como um exemplo de análise discriminante linear”. Às vezes, é chamado de conjunto de dados da íris de Anderson porque Edgar Anderson coletou os dados para quantificar a variação morfológica das flores da íris de três espécies relacionadas. Duas das três espécies foram coletadas na Península de Gaspé, “todas do mesmo campo, colhidas no mesmo dia e medidas ao mesmo tempo pela mesma pessoa com a mesma aparelho”. O conjunto de dados consiste em 50 amostras de cada uma das três espécies de Iris ( Iris setosa, Iris virginica e Iris versicolor). Quatro variáveis foram medidas em cada amostra: o comprimento e a largura das sépalas e pétalas, em centímetros. Com base na combinação dessas quatro características, Fisher desenvolveu um modelo discriminante linear para distinguir as espécies umas das outras (Fonte: https://pt.wikipedia.org/wiki/Conjunto_de_dados_flor_Iris).
#Pacotes utilizados
library(klaR)
library(psych)
library(MASS)
library(ggord)
library(devtools)
# Obtendo os dados
data("iris")
str(iris)
## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
#Criando um gráfico de dispersão para as quatro primeiras variáveis numéricas.
pairs.panels(iris[1:4],
gap = 0,
bg = c("red", "green", "blue")[iris$Species],
pch = 21)
Se pode ver no gráfico, diagrama de dispersão, histograma e valores de correlação. Em seguida, se quer criar os melhores grupos de separação com base nessas espécies. Vai ser criado um conjunto de dados de treinamento (training) e um conjunto de dados de teste (test) para fins de previsão e teste. 60% do conjunto de dados usado para fins de treinamento e 40% usados para fins de teste.
#Separação dos grupos
set.seed(123456789)
ind <- sample(2, nrow(iris),
replace = TRUE,
prob = c(0.6, 0.4))
training <- iris[ind==1,]
testing <- iris[ind==2,]
# Análise discriminante linear
linear <- lda(Species~., data=training)
linear
## Call:
## lda(Species ~ ., data = training)
##
## Prior probabilities of groups:
## setosa versicolor virginica
## 0.3043478 0.3478261 0.3478261
##
## Group means:
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa 5.035714 3.460714 1.475000 0.2357143
## versicolor 5.896875 2.731250 4.215625 1.3125000
## virginica 6.553125 3.034375 5.509375 2.0281250
##
## Coefficients of linear discriminants:
## LD1 LD2
## Sepal.Length 0.5851444 0.4528399
## Sepal.Width 2.0168111 -2.5347510
## Petal.Length -1.8643109 0.6817960
## Petal.Width -3.2514946 -2.7190523
##
## Proportion of trace:
## LD1 LD2
## 0.9855 0.0145
attributes(linear)
## $names
## [1] "prior" "counts" "means" "scaling" "lev" "svd" "N"
## [8] "call" "terms" "xlevels"
##
## $class
## [1] "lda"
# Histograma empilhado para valores de funções discriminantes.
p <- predict(linear, training)
ldahist(data = p$x[,1], g = training$Species)
Com base no conjunto de dados de treinamento, 30,4% pertencem ao grupo setosa, 34,7% pertencem aos grupos versicolor e 34,7% pertencem aos grupos virginica. A primeira função discriminante é uma combinação linear das quatro variáveis. A porcentagem de separações alcançada pela primeira função discriminante é 98,55% e a segunda é 1,45%.
Esses histogramas são baseados na primeira função discriminante. É evidente que não há sobreposições entre primeira e segunda e primeira e terceira espécies. Contudo, alguma sobreposição é observada entre a segunda e a terceira espécie.
ldahist(data = p$x[,2], g = training$Species)
Os histogramas baseados na segunda função discriminante são sobrepostos, não bons.
ggord(linear, training$Species, ylim = c(-10, 10))
Os gráficos Biplot acima são baseados nas duas funções discriminantes. Setosa se separou muito claramente e alguma sobreposição foi observada entre Versicolor e Virginica.
Com base nas setas, a largura e o comprimento da sépala explicaram mais para Setosa, a largura e o comprimento da pétala explicaram mais para Versicolor e Virginica.
Finalmente, é possível ver as classificações corretas e classificações de erros.
p1 <- predict(linear, training)$class
tab <- table(Predicted = p1, Actual = training$Species)
tab
## Actual
## Predicted setosa versicolor virginica
## setosa 28 0 0
## versicolor 0 30 1
## virginica 0 2 31
sum(diag(tab))/sum(tab)
## [1] 0.9673913
p2 <- predict(linear, testing)$class
tab1 <- table(Predicted = p2, Actual = testing$Species)
tab1
## Actual
## Predicted setosa versicolor virginica
## setosa 22 0 0
## versicolor 0 18 0
## virginica 0 0 18
sum(diag(tab1))/sum(tab1)
## [1] 1
No conjunto de dados de treinamento, a classificação correta total é 30+30+29=89. A acurácia do modelo é aproximadamente 100%.
No conjunto de dados de teste, a classificação correta total é 20+19+20=59. A acurácia do modelo é aproximadamente 96,72%.
Assim, o Histograma e o Biplot fornecem informações úteis para interpretações e, se não houver uma grande diferença nas matrizes de covariância do grupo, a análise discriminante linear funcionará muito bem. A Análise Discriminante Linear não é útil para resolver problemas não lineares.